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White's density matrix renormalization group (DMRG) method has been ap- 
plied to the one-dimensional Ising model in a transverse field (ITF), in order 
to study the accuracy of the numerical algorithm. Due to the exact solubil- 
ity of the ITF for any finite chain length, the errors introduced by the basis 
truncation procedure could have been directly analysed. By computing dif- 
ferent properties, like the energies of the low-lying levels or the ground-state 
one- and two-point correlation functions, we obtained a detailed picture, how 
these errors behave as functions of the various model and algorithm param- 
eters. Our experience with the ITF contributes to a better understanding of 
the DMRG method, and may facilitate its optimization in other applications. 
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I. INTRODUCTION 



In the last three years we have witnessed a breakthrough in the numerical analysis of 
one-dimensional (ID) quantum lattice models. This considerable progress was due to the 
invention of the density matrix renormalization group (DMRG) method by S. White.0 Recent 
applications of the method has demonstrated its extreme efficacity and versatility. By now 
the DMRG has become one of the leading numerical tools in the study of most ID quantum 
spin and electron problems of current interest. 

There are several numerical methods to obtain the low-energy states of a given quantum 
Hamiltonian. Exact diagonalization algorithms are able to compute the ground state and 
the lowest excited states with a precision of more than 12 digits. The attainable system 
sizes, however, are rather limited due to memory restrictions. On the other hand, stochastic 
methods, like the diverse variants of the quantum Monte Carlo method, are capable of 
treating systems with hundreds of sites - at the price of reducing the precision of the obtained 
results. In this case, it is rather the computation time that limits the applications. Use 
of renormalization group procedures, such as the DMRG, could be the way out of this 
dilemma. The key idea is to gradually increase the system size, and, at the same time, 
systematically truncate the Hilbert space by keeping only those degrees of freedom that 
are really "important" for an accurate representation of the desired state. The main point, 
and this is where the DMRG differs drastically from the preceding numerical RG methods, 
is how to choose the most important degrees of freedom to minimize the error caused by 
discarding the other "unimportant" ones. 

The DMRG is an iterative algorithm to build up the lattice to the desired length and 
find approximants to the target state (the ground state or an excited state), using only a 
limited number of basis states. The total system is divided into two parts, the block and the 
environment. In each step of the procedure the block is increased by adding one lattice site, 
but only the most probable states, obtained from the reduced density matrix of the block, 
are kept. The environment is composed from the block of the previous step, and its role is 
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to embed the block into a larger system, when the density matrix is formed. Application 
of the reduced density matrix in the truncation process, and not simply keeping the lowest 
energy states of the block as in previous RG techniques^ is the key ingredient of the DMRG 
method, since, as was shown by White,0this is the way to minimize the error introduced into 
the representation of the target state. Although, during the algorithm the length of the total 
system increases gradually, the dimension of the Hilbert space is always kept manageable 
by the truncation process, and systems with relatively large size can be studied. 

The DMRG has been successfully applied to various ID and coupled chain problems, 



problems^ or the two-chain Heisenberg and Hubbard models.Q Promising implementations of 
the algorithm to compute dynamical properties! of ID systems or to simulate 2D latticesi§ 
have also been reported. At least for quasi ID problems, its advantage over the standard 
numerical procedures has become evident. While the available system sizes are comparable 
to those of the Monte Carlo method, the precision of the computed quantities seems better 
by several orders of magnitude. 

The DMRG works especially well when the system is subject to open boundary condition. 
In the case of periodic boundary condition, on the other hand, errors are definitely stronger. 
Moreover, the conservation of momentum cannot be directly built into the method, so re- 
striction of the computation to a certain momentum sector is not possible. These drawbacks 
has led to the fact that a large part of current applications of the method treat open systems. 

Although this is a renormalization group procedure, contrary to the naive expectations, 
results seem more accurate if the system is away from criticality, i.e., when the model 
possesses a finite spectral gap in the thermodynamic limit. Other parameters of the model 
in study, like the number of degrees of freedom at a single site or the range of interactions 
also have drastic effects on the precision. 

For any numerical method, it is of essential importance to understand, how errors evolve, 
when parameters of the procedure or those of the model in study are varied. So far, the 
accuracy of the DMRG method has only been tested on small clusters where numerically 
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exact results were known from exact diagonalization, or alternatively, large lattice DMRG 
data were extrapolated to the infinite chain limit, and these extrapolated values have been 
compared to rigorous results, available only in the thermodynamic limit. In many of the 
applications, however, at most the small lattice exact data are available for testing. Nev- 
ertheless, large lattices, even with hundreds of sites, are used to draw conclusions on the 
behavior of the system. In many cases, it is really a difficult problem to reliably estimate 
the precision as the system size increases. 

Our main goal in this paper is to study the general trends of numerical errors during the 
DMRG algorithm, especially when the length of the system is in the range of the typical 
applications, i.e., well beyond the limit where exact diagonalization data are attainable. 
For this purpose, we apply the method to an exactly solvable system, the Ising model in a 
transverse field (ITF). The main advantage of the ITF model as a test system is that its 
energy spectrum and ground-state correlation functions can be calculated by simple tools 
for any finite chain length, even when the chain is subject to open boundary condition, as 
in the standard DMRG applications. The ITF is also a model with a second order quantum 
phase transition, so the effect of the appearance of criticality on the accuracy can also be 
analysed. 

We carried out a detailed DMRG study of the ITF, varying several parameters of the 
model and the numerical procedure, like the strength of the magnetic field, the chain length, 
the number of states kept, the number of iterations and the number of target states. Beside 
the energies of the low-lying states, we also computed different one- and two-point correlation 
functions. 

The setup of the paper is as follows. In Sec. II we briefly summarize the exact solution of 
the ITF with open boundary condition. Section III is devoted to the details of the DMRG 
algorithm, while Sec. IV contains the analysis of the observed trends in the numerical errors. 
Section V is a summary of our main conclusions. 
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II. THE ID ITF WITH OPEN BOUNDARY CONDITION 



The Ising model in a transverse field (ITF) on a one- dimensional chain of N sites is 
defined by the following Hamiltonian 

N-l N 

w = -E« + i-7E<. (i) 

n=l n=l 

where a a , a = x,y,z are the Pauli matrices, 7 is the transverse magnetic field applied 
in the ^-direction, and the chain is subject to open boundary condition, as in our DMRG 
implement at ion . 

In the thermodynamic limit, the ITF possesses a second order phase transition that 
takes place at 7 = 1, where the energy gap vanishes as N — > 00. For 7 < 1 (ordered phase), 
the ground state is doubly degenerate and there is a long-range order: the ground-state 
correlation function pf = (c^a^ +l ) tends to a finite value as I — > 00. On the other hand, for 
7 > 1 (disordered phase), the ground state is unique and pf — > as I — > 00. The gap is 
finite on both sides of the transition point. 

The ITF in Eq. (|l|) is exactly solvable for any chain length N. Details of the calculation 
can be found, e.g., in Ref. |Tl|. Here we only summarize the main steps and the necessary 
formulae. Introducing Fermi operators via the Jordan- Wigner transformation 

c n = exp (^ri ^ tfvjj (2) 
the Hamiltonian 7i reduces to a quadratic form 

N N-l 



H = -2 7 £ 4c n - (4c„ +1 + 44 +1 + h.c.) + iV 7 . (3) 

n=l n=l 

This can be diagonalized directly by a Bogoliubov transformational 

Vk = £ + , (4) 

leading to the final diagonal form 

n= E m(vlv k -l), (5) 



k allowed 
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A(fc) = 2 (l +7 2 + 2 7 cosA;) . (6) 
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The TV-element real vector <p k is a solution of a set of linear equations,0 and turns out to 
be 

6 k = An (sin kn - cos kn | , (7) 

\ 7 + cos k J w 

where n = 1, 2, . . . , N and is a normalization constant. i\) k can be expressed from (p k as 

^ = A' N ( 7 0* + <f> k n+1 ) , (8) 

where A' N is another normalization constant, and formally ip% = 0. 

Reflecting the fact that the chain is subject to open boundary condition, the allowed k 
modes in the summation of Eq. (|6]) cannot be written in a simple form for general 7, but 
are determined through a trigonometric equation 

sin[fc(iV + 1)] = _ 1 

sin(A;iV) 7' 1 ' 

In any case, the total number of independent modes is equal to the number of sites N. By 
convention, the allowed k values are the roots k = k + iki, whose real part ko is in the 
interval < k < it. When 7 > 7 c (A r ) = N/(N + 1) all the roots are real. For 7 < 7 C (JV), 
however, a complex (localized) solution becomes possible with ko = ir and k% > 0. Note that 
7 c (iV) — > 7 C = 1, the phase transition point, as iV — >• 00. The imaginary part k\ can be 
obtained by solving the equation 

sinh[fci(iV + 1)] _ 1 



smh(k\N) 7' 

and the corresponding <f> vector is0 



(10) 



<tC = A^-lJ 1 - 1 sinh [(TV — n + , (11) 

where, again, n = 1,2, ... ,N and A^ is a normalization constant. ip l ° c is obtained again 
through Eq. (|8|). The elements of (f) l ° c and ip l ° c are exponentially small far from the left 
and right chain ends, respectively, which indicates, by Eqs. (f|), that the complex mode is 
localized to the chain ends. The energy associated with this localized mode is 
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A(7r + ifci) = 2(l+7 2 - 27coshfci) . (12) 



When N is large, Eq. (|TU|) is readily solved to yield k\ ~ In I/7. By substitution, this gives 
A(-7r + iki) ~ 0, i.e., it is a mode of zero energy, leading to the double degeneracy of the 
ground state sector when 7 < 1. 

The ground state energy of the model, the fermionic vacuum, is expressed as 

E, = \ £ A (k), (13) 

k allowed 

that can be written in a closed form only at the transition point 7 = 1, where it reduces^ 
to 

E " = 1 " cosec (4FT2) ■ (14) 

The first excited state of the model is always determined by the wave number k, whose 
real part kn is closest to 7r; for 7 > 7 C it is a real root and the energy gap converges to 
2(7 — 1) as N — > 00. For 7 < j c , this is the localized mode, and the first excited state 
becomes asymptotically degenerate with the ground state. It is the second excited state 
that constitute the real (finite) energy gap in this case. 

Finally, ground-state one- and two-point correlation functions can also be calculated, by 
the method described in details by Lieb et aZ.tl. The only interesting one-point function is 
M. z n = (cr^), which can be expressed as 

M z n = l -G nni (15) 

where G nm is defined by 

G nm = - £ V£#r (16) 

k allowed 

Note that even in the ordered state }A x n = (a*) = 0, since the Hamiltonian, and thus 
any finite lattice ground state, is invariant to the transformation a* — > —o x n . The same 
holds for cr^, so that M. y n = (a^) = 0, too. As for the two-point functions, p x nm = (cr^a^), 
Pnm = ( a n a m) an d Pnm = ( a n a m ) i one arrives to the formulae 



P X nm = 1/4 



Pnm 

= 1/4 



G n ,n+1 Gn,n+2 



G m —l,n+l 



Gm,n 



Gr. 



G m —l,m 
G>i+l,m— 1 



G m ,m—1 



(17) 



(18) 



and 



Pr, 



1/4 (G nn G mm G mn G nm j , 



(19) 



where n < m. 

In order to be able to compare the approximate results, produced by the DMRG method, 
with the exact ones for a certain value of 7 and chain length N, the exact quantities were 
calculated by an independent numerical process, using either double precision Fortran rou- 
tines or the Mathematica software that allows computations with arbitrary precision. We 
solved the nonlinear equations (|]) and fllQ|) to obtain the allowed set of wavenumbers, and 
used these values to express the associated 4> k and i\) k vectors. Then the ground-state and 
excited state energies and the correlation functions were calculated numerically according 
to the above formulae. All the results obtained in this way were precise to at least 14 digits. 



III. THE DENSITY-MATRIX RENORMALIZATION GROUP METHOD 

White's density matrix renormalization group method (DMRG) is a numerical real-space 
renormalization group procedure, in which the effective size of the system increases gradually, 
while the dimension of the associated Hilbert space remains constant, due to a systematic 
truncation process. Since the method is well described in the original papers,^ we only 
present a brief summary here. 

The basic object of the method is the system block Bi that consists of I lattice sites. All 
the necessary operators A4 of this block (e.g., cr", a = x,y,z and 1 < n < I) are stored as 
M x M matrices. Note that the real dimension of a block of I sites is d l , where d is the 
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number of states at a single site, and M « d , so the representation of the block matrices 
is only approximative. In each step of the algorithm, a new single site • is added to the 
existing block Bi, and operators of the resulting system B\» of dimension Md are formed as 
tensor products from the matrix representations of the corresponding operators of the two 
constituting parts. Bi* is then renormalized, first by carrying out an appropriate unitary 
transformation, then by truncating its degrees of freedom from Md to M. Only the most 
important M states are kept, and the remaining M(d — 1) ones are discarded. The resulting 
system of dimension M and length I + 1, denoted by -B;+i, is then used iteratively in the 
subsequent step of the algorithm. Choosing the states kept, i.e., finding the appropriate 
unitary transformation, is the most crucial point of the renormalization group procedure. 

Previous applications of the renormalization group technique, where the M lowest energy 
states of Bi» were kept, led to disappointing results for many model systems. As was 
observed by several authors,0 this was due the interaction having been neglected during 
the renormalization process between Bi* and the rest of the larger (in many cases infinite) 
system of what it makes a part. Bi* was renormalized in a way that an unnatural open 
boundary condition was forced on both of its ends. 

In order to avoid this problem, White's method embeds Bi* into a superblock, and uses 
the eigenstates of this larger system to carry out the renormalization of Bi*. Choosing 
the superblock configuration as Bi • •-B/ 2 , where Bf is the reflection of Bi, its effective 
Hamiltonian, built up from the matrix representations of the necessary operators of the two 
blocks and the two single sites, is diagonalized to obtain the desired (target) state Even 
though, the target state can be a linear combination of more states, by targeting only one 
state, the renormalized block states are more specialized for representing that single one, 
and fewer of them are needed for a given accuracy. 

Having the target state expressed on the superblock as 

Md 

*= E *yl»*>|j>, (20) 

where i and j label the Md states of Bi» and its surroundings •-B/ J , respectively, the reduced 
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density matrix of the subsystem Bi» is formed as 



Pii' ^ t ^ ij^ i'j ■ 



(21) 



j 



As was shown by Whitejil the error introduced into the representation of the target state 
by the truncation to M states is minimized, if the new basis that one changes to is the 
basis that diagonalizes the density matrix pu>. The renormalization of the tensor product 
operators of 



is thus carried out by the M x Md transformation matrix O, whose rows are composed from 
the M density matrix eigenvectors, associated to the M largest eigenvalues cu a , a = 1, • • • , M. 

The simplest version of the DMRG method, the infinite lattice algorithm, starts with 
four lattice sites, i.e., from the superblock configuration Bi • •Bf. In each step, the total 
length of the chain increases by 2. Measurements of the interested quantities are made 
after the calculation of the target state in each step, and the whole process is continued 
until the results converge satisfactorily. This method is especially suitable to yield, with 
minimal computational efforts, rather precise estimates of bulk properties, like the ground- 
state energy density or ground-state correlation functions. 

The usual measure of the accuracy of the truncation to M states is the deviation of 
the sum of the density matrix eigenvalues associated with the states kept Pm = Ha=i w « 
from unity. Clearly, in the extreme case when the discarded states have zero weight, i.e., 
— = ■ ■ ■ = ujMd = so that Pm = 1, they are not required to represent the 

target state and no error has been committed. There is, however, another source of 
error in the DMRG procedure, namely that the B\» subsystem is only embedded into an 
approximate superblock when it is renormalized, and not into the, a priori unknown, exact 
environment. The two kinds of errors, the "truncation" error and the "environment" error 
are not simply additive. The latter one, however, can be reduced in an iterative manner, 
using the so called finite lattice algorithm. 



At -> OAiO\ 



(22) 
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The finite lattice algorithm starts by building up the lattice to the desired length N, 
using the infinite lattice algorithm. Then the superblock configuration is modified to B\ • 
•B§_ 2 _i, so that the total length remains always N, and the blocks Bi, I — 1, • • • , N — 3, are 
recomputed. The process when all B^s are recomputed constitute the cycle of the algorithm. 
Since at each step from this time on, the block to be renormalized is part of a system of the 
desired length N, the environment becomes more precise, and a considerable improvement 
of the results, corresponding to the TV-site system, is achieved. These results can be then 
used to carry out a systematic finite-size scaling analysis, or to study, e.g., boundary effects 
in the finite chain. 

IV. NUMERICAL RESULTS 

In order to determine the accuracy of DMRG method, numerical calculations on the 
ITF, using both the infinite and finite lattice algorithms, were performed. Errors in various 
quantities, such as the ground-state and first excited-state energies Eqs and -Eixs, resp., 
the one-point correlation functions K\ z n and M. x w and the two-point correlation functions 
p\ and pf were monitored. In the case of the infinite system algorithm, our main concern 
was how the errors of the energies depend on the system size N, when other parameters are 
kept constant. In the finite lattice algorithm, on the other hand, the length of the chain 
was kept fixed, and the effect of introducing additional cycles in the process was analyzed. 
When excited states were computed, the number of target states was also varied. 

Since the ITF possesses a second order phase transition at 7 C = 1, it was expected that 
many of the above quantities have different behavior, whether 7 is equal to the critical value, 
or greater or less than 7 C . Therefore, we investigated the three different regions of the phase 
diagram by choosing three different values for 7, namely, 7 = 0.5 < 7c , 7 = 1 = j c , and 
7 = 2 > 7c . 
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A. Energies 



Let us first consider the error of the ground-state (GS) and first excited-state (1XS) 
energies. Using the infinite lattice algorithm, we built up a chain of length N = 300 and 
kept states up to M = 48. Except one example which will be discussed below, we found that 
targeting one state alone is by orders of magnitude more precise that targeting several states 
together. Hence, unless stated otherwise, the results to be presented here were obtained in 
a way that the ground state and the first excited state were targeted separately. 

The errors E'(dmrg) — E (exact) were always found to be positive, satisfying the statement 
that the DMRG is a variational method, which gives (at least for the GS energy) an upper 
bound estimate.! The relative errors 5E = [£(dmrg) — E (exact)]/ |J5(exact)| as a function 
of iV for various values of M and 7 are shown on log-log scale in Fig. [l|-|3|. 

Until the algorithm keeps all states, i.e., for iV < 2(lnM + 1), the DMRG is exact and 
only the machine's numerical inaccuracy is seen, limiting the precision on the order of 10~ 14 . 
For longer chains, however, the errors associated with the reduction of degrees of freedom 
also come in. 

For the off-critical values, 7 = 2 and 0.5, the ground state error 5Eq§ shows practically no 
size dependence. It depends, however, crucially on the value of M. The behavior is the most 
clear for 7 = 2, where the infinite system ground state is unique with a gap above it. Even 
a small value of M, M ~ 10, is enough to reach the machine's precision limit [Fig. |(a)]. 
There is, on the other hand, an interesting size dependence in the error of the first excited 
state SEixs, especially for greater M values [Fig. 3(b)]. A maximum evolves, above which 
SEixs begins to dwindle again. This can be understood, however, since the approximate 
excited state, yielded by the DMRG process, is not necessarily exactly orthogonal to the real 
ground state, and can have a finite overlap with it. This fact leads to an overcompensation 
of the otherwise increasing positive absolute error for long chains. 

For 7 = 0.5, where the ground state is asymptotically doubly degenerate, the behavior 
is more subtle [Fig. 0(a,b)]. The splitting of the two lowest levels decreases exponentially 
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as N increases. When this difference goes below the machine's precision limit ~10 -14 , 
the algorithm is incapable to further resolve the two levels, and the target state that it 
yields is a linear combination of the two exact eigenvectors with random coefficients. As a 
consequence, in each step of the DMRG algorithm the target state changes unpredictably, 
leading to the loss of optimality of the truncation process, and hence, to considerably higher, 
rather scattered error rates. As Fig. |2|(b) shows, the DMRG can even " lose" the target state 
above a certain chain length. In a try to avoid these problems, we also computed E\xs 
by choosing the target state to be a linear combination of the ground-state and the first 
excited-state vectors. This proved to be an improvement in the range where otherwise the 
excited state was lost; for smaller values of N, however, the errors were considerably larger 
[Fig. |(b)]. 

At the critical point 7 = 1, the errors were found to be by several orders of magnitude 
larger than at the off-critical values of 7. This is in accordance with the findings of Ref. 



15| : the correlation length of the model is one of the most significant factors that influence 
the precision of the DMRG method. Moreover, the curves in Fig. |3|(a,b) have a clear size 
dependence. When M is kept fixed, the errors can become larger by 4 to 6 orders of 
magnitude, as the chain length approaches N = 300. For smaller M, the clear downward 
curvature seen in the log-log plot indicates that the errors converge to a finite value at large 
lengths. For larger M values, however, this convergence is much slower, and the analysis is 
made more difficult by the appearance of a crossover effect, which change the behavior in 
the small N region, and makes the curves more flat there. While for M = 16 the crossover 
size is around iV ~ 10 (and hence unobservable) , for M = 32 it is at iV ~ 100, and for 
M = 48 it is at iV ~ 250, showing that the crossover size scales for larger lengths as more 
and more states are kept. 

A possible interpretation of this crossover effect can be obtained by recalling that there 
are two different sources of errors in the DMRG method (see Sec. III). For small M, this is 
clearly the "truncation" error that dominates. Curves with M = 4, 8, 16 in Fig. |3](a,b) basi- 
cally show the size dependence of this type of error alone. The effect of the "environment" 
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error only shows up for large enough M and small enough N values, when the truncation 
error is strongly reduced. The environment error approaches its saturation earlier than the 
truncation error, as is seen from the M = 48 curves, where the environment error dominates. 
However, since the truncation error increases several orders of magnitude in the small N 
regime, it always exceeds the environment error for large enough sizes. This produces a 
rather sharp break in the curves in the log-log plot at the crossover size, such as seen for 
M = 32 in the figure. 

For comparison, and to reduce the environment error, computations using the finite 
lattice algorithm were also performed. Fixing the chain length at iV = 100, the iteration 
process was repeated until the desired energies converged. For 7 = 2, the relative errors are 
plotted in Fig. f|(a,b) as a function of the sum of the discarded density matrix eigenvalues 
1 — P m . While there is practically no improvement in the ground-state energy, the first 
excited-state energy becomes much more precise, when further cycles are carried out, and 
full convergence is reached only after the 1 = 3 iteration. At the critical point 7 = 1, SEqs 
and SEixs behave similarly: two cycles are needed to get rid of the environment error [Fig. 
§(a,b)]. 

It is seen in the figures that the 1=1 cycle data (the infinite lattice algorithm results) 
do not fit onto a straight line on the log- log plot. Points from the fully converged cycles, 
however, do so nicely. The slope of the fitted line was found to be very close to unity in all 
cases, indicating that the error is proportional to the discarded weights, i.e., 

5E = const -(1-P m ), (23) 

where the constant can depend on the model parameters and the system size. We emphasize, 
however, that this form only holds for the converged energies. Extrapolating the infinite 
lattice algorithm (J = 1) data by this formula to the 1 — P m — * (M — ► oo) limit can yield 
false results. 
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B. One-point correlation functions 



Expectation values of local operators in the ground state were computed using the finite 
lattice algorithm at a fixed chain length N = 100. Both M.^ and M. z n (n=l, . . . ,100) were 
measured, and the errors 5A4% = 7V1"(dmrg) — .M" (exact), a = x, z were calculated. Note 
that for Ai^, the exact values are zero, as it was detailed in Sec. II. 

At the critical point 7 = 1, our results are presented in Fig. |](a,b). The 1=1 cycle 
(infinite lattice algorithm) produces an error which depends significantly on the position in 
the chain n. Results are more precise around the middle of the system. Note that spatial 
variations can reach several orders of magnitude for larger M values, as e.g. for the 1=1, 
M = 16 curve in Fig. ||(a). There is also a sudden improvement in accuracy very near to 
the chain ends, but this is believed to be an anomaly of the ITF model and not a general 
feature of the DMRG technique. (The exact one-point functions of the model, when open 
boundary condition is used, show a strong boundary effect close to the ends, and this seems 
to influence the errors too.) 

In the case of Ai z n [Fig. ^(a)], additional iteration cycles (I > 2) considerably improve 
the precision by decreasing the environment error. Although some fluctuations may still be 
present (especially for large M) in the 1 = 2 data, the error becomes more or less constant 
and spatially more homogeneous for I > 2. The convergence found with respect to the 
number of cycles is similar to that of the energies. 

For .M* [Fig. P(b)], where the exact values of the magnetization are zero, the situation 
is different. Additional cycles, and increasing the value of M, rather unexpectedly, make 
the data less precise. The sign of the error also changes: while for the 1 = 1 cycle SAi^ is 
negative, for / = 3 it always turns out to be positive. (Note that in the figure the absolute 
values of the errors are plotted.) The rapid oscillations seen in the / = 2, M = 8 curve are 
due to the fact that the error fluctuates between positive and negative values. We believe 
that for this pathological case the observed errors stem from the numerical inaccuracy of the 
diagonalization subroutines, and not from the algorithm itself. Test runs on short systems 
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keeping all states in the blocks, i.e., when the DMRG is numerically exact, produced the 
same qualitative picture of 

It might have been thought that monitoring the magnitude of errors of quantities, whose 
exact results are a priori known, like M. x n in our case, could yield information on the precision 
on other, a priori not known quantities, like M. z n - We see, however, that the errors of the 
two one-point functions behave completely differently without any obvious correlation, so 
the knowledge of the accuracy of one of them cannot be used to draw predictions on the 
other one. 

For the 7 = 2 case (figure not presented), all the iteration cycles gives the same result 
for M. z n . There is practically no n-dependence, and the algorithm reaches the border of 
numerical inaccuracy very soon, in accordance with what was found for the ground-state 
energy. For M = 8, the relative errors scatter on the scale of 10~ 9 , and this precision could 
not be improved by increasing M. The behavior was found to be very similar in the ordered 
phase, at 7 = .5, the only difference is that more cycles were needed for the full convergence. 

C. Two point correlation functions 

Measurements of the two-point correlation functions were carried out similarly to the one- 
point correlation functions, at a fixed length N = 100. Following White's recipe, pf = p® n+l , 
a = x,z, was measured so that the points n and n + 1 were positioned symmetrically to the 
middle of the chain.0 This special allocation of the points assures that, at least for the short- 
range correlations, end effects are strongly reduced. Although, in most applications the bulk 
correlation functions are of interest, since our aim was to test the DMRG algorithm itself, 
we compared the numerical data with the exact finite-lattice results obtained at iV = 100. 

At the critical point, 7 = 1, the absolute value of the error of the correlation functions 
\Spf\ = \pf (dmrg) — pf (exact)| is plotted in Fig. 0(a,b) for a = z and x. For small M values, 
the curves are rather smooth with a moderate dependence on I. Convergence with respect 
to the number of cycles I is reached when the ground-state energy converges. For larger 
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values of M, the errors seem to fluctuate rather irregularly, especially for \Spf\, for which 
several cusps evolve, where the error rate drops abruptly. This strange feature, however, is 
once again an artifact of plotting the absolute value of the error on a logarithmic scale. For 
each I where the cusp appears, we found that Spf changes sign, similarly to the behavior of 
the corresponding }A x n curve, and this sign fluctuation causes the strange-looking shape of 
the curves. 

Disregarding the fact that the sign of the errors is not fixed, there is a clear tendency 
that in the 1 = 1 cycle (infinite lattice algorithm) the long-range correlators are less precise. 
This is in complete accordance with the behavior of the one-point functions: operators are 
represented less accurately moving towards the chain ends because of the environment error. 
More cycles improve the situation, especially for 8pf [Fig. 0(a)], where the errors for large 
/ decrease by 3 orders of magnitude and, quite unexpectedly, the long-range correlators 
become more precise than the short-range ones. On the other hand, there is no similar 
change in the tendency of the curves for Spf [Fig. 0(b)]. Although the improvement is the 
most significant for large I, long-range correlators remain less accurate. 

For the noncritical values 7 = 2 and .5 (figure not presented), Spf and Spf show practically 
no /-dependence. There is only a moderate change in accuracy for very small and very large 
values of I, but this does not exceed an order of magnitude. The average accuracy of the 
correlation functions is considerably worse than that of the energies. The limit of the relative 
precision we could achieve by increasing M was not better that 10 -9 , similarly to what was 
found in the case of the one-point functions. The effect of carrying out more cycles was also 
similar, the correlation functions converged exactly when the energies reached their limit 
values. 

V. SUMMARY 

In the present paper, we have analysed in detail the accuracy of White's density matrix 
renormalization group method, by applying it to the one- dimensional Ising model in a trans- 
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verse field. Due to the exact solvability of this model, the exact and the numerical results 
could have been directly compared. We varied several parameters, either in the numerical 
algorithm or in the model, and obtained a rather detailed picture how the accuracy of the 
DMRG approximation depends on these parameters. 
Our main results are summarized as follows: 

(i) The DMRG yields an extremely precise value for the ground-state energy, especially 
when the model is far from the critical point, and the ground state is unique. Excited state 
energies or the ground state energy of a critical system can be obtained with much less 
accuracy. 

(ii) Targeting exclusively one of the excited states is, in general, more precise than 
targeting it together with the ground states or other low-energy states. This, however, can 
become unstable when the levels are too close to each other, like in the case of asymptotic 
degeneration of the ground state. When this is expected to happen, more states must be 
targeted together. 

(hi) Carrying out only the 1 = 1 iteration cycle (infinite lattice algorithm) leads to a 
significant "environment" error. Its effect is the most pronounced, when lots of states are 
kept, i.e., when M is large, so the "truncation" error is relatively small. The two types of 
error produce a crossover effect as N increases in the critical case. The environment error 
dominates in the small N large M regime. 

(iv) Although the chain length dependence of the errors is strong in the critical case, both 
the truncation and the environment errors seem to converge to finite values as N increases. 
In the practical range of applications 10 < N < 300, however, the errors, especially the 
truncation error, become larger by 2-3 orders of magnitude. 

(v) For large M values, the finite lattice algorithm (7 > 2) improves the results consid- 
erably by decreasing the environment error. This is seen not only in the energies, but in the 
one- and two-point correlation functions. When the data are fully converged (but not after 
the first 1=1 cycle), the errors are nicely proportional to the sum of the discarded density 
matrix eigenvalues 1 — P m . This fact can be used to make an extrapolation to the M — > oo 
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limit. 

(vi) Accuracy of the correlation functions is always worse than that of the energies. In 
the 1 = 1 cycle, the error in the representation of the local operators become larger as 
one moves outward from the chain center. When two-point functions are computed in the 
usual way, i.e., symmetrically to the chain center, this leads to the fact that the long-range 
correlators are less precise. Additional iteration cycles make the errors smaller and their 
spatial dependence more homogeneous. 

Although the one-dimensional ITF is a rather simple many body system, we believe that 
most of our findings hold equally well for other, more complex quasi-one-dimensional lattice 
problems. We hope that the above results contribute to a better understanding of the DMRG 
procedure, and provide a direct help in optimizing the algorithm in other applications. 
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FIGURES 

FIG. 1. Relative errors of the (a) ground-state energy, and (b) first excited-state energy as a 
function of the chain length N for different values of M. Data obtained by the infinite lattice 
algorithm (I = 1) in the disordered regime at 7 = 2. 

FIG. 2. Same as Fig. |], but in the ordered regime at 7 = 0.5. Curves labelled by TS = 2 show 
the case, when the excited state was targeted together with the ground state. 

FIG. 3. Same as Fig. [j], but at the critical point 7 = 1. Curves with M = 4,8, 16 are dominated 
by the truncation error, while those with M = 48 by the environment error. The M = 32 curves 
show a crossover between the two types. 

FIG. 4. Relative errors of (a) the ground-state energy, and (b) first excited-state energy as 
a function of the sum of the discarded density matrix eigenvalues 1 — P m for different values of 
M (shown as labels). Data obtained by the finite lattice algorithm, carrying out the I = 1,2,3 
iteration cycles, and the chain length is fixed at N = 100. Disordered regime at 7 = 2. 

FIG. 5. Same as Fig. but at the critical point 7 = 1. 

FIG. 6. Errors of the one-point functions (a) Ai^, and (b) M% as a function of the position 
in the chain n. Data obtained by the finite lattice algorithm, carrying out the I = 1,2,3 iteration 
cycles, and keeping M = 8 or 16 states. The chain length is fixed at N = 100. Critical point 7 = 1. 

FIG. 7. Errors of the two-point functions (a) pf , and (b) pf as a function of I. Data obtained 
by the finite lattice algorithm, carrying out the I = 1, 2, 3 iteration cycles, and keeping M = 8 or 
16 states. The chain length is fixed at N = 100. Critical point 7 = 1. 
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